Intro to BioEnergeticFoodWebs

by Chris Griffiths, Eva Delmas and Andrew Beckerman, Dec. 2020.

7.7 μs

This document follows on from 'Getting started', 'Basic Julia commands' and 'Differential Equations in Julia' and assumes that you're still working in your active project.

This document introduces the BioEnergeticFoodWebs.jl and EcologicalNetworks.jl packages. It demonstrates how to run the BioEnergetic Food Web (BEFW) model, how to vary variables of interest (e.g., productivity) and construct experiments designed to investigate the effect of different variables on population and community dynamics.

For those that are unfamilar with the BEFW and it's application in Julia, we advise checking out the MEE paper before we start. Remember, the BEFW model is also based on a system of differential equations and is solved using the same engine as the DifferentialEquations.jl package.

7.0 μs

Load packages

3.0 μs

You'll need the following packages for this tutorial:

2.7 μs
36.4 s

The JLD2.jl package will be useful later as it allows you to directly export and load a BEFW output object. Let's also set a random seed for reproducibility:

3.1 μs
122 ms

Preamble

One of main advantages of running food web models in Julia is that simulations are fast and can be readily stored in your active project. With this in mind, make a new folder in your project called out_objects (right click > New Folder). Alternatively, you can create an out_objects folder directly using mkdir().

4.9 μs
"example_folder"
101 μs

Running the BEFW

There are four major steps when running the BioEnergetic Food Web model in Julia:

  1. Generate an initial network

  2. Fix parameters

  3. Simulate

  4. Explore output and plot

Initial network

Before running the BEFW model, we have to construct an initial random network using the niche model. The network is characterised by the number of species in the network and its connectance value. Here, we generate a network of 20 species with a connectance value of 0.15:

11.6 μs
20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0
36.2 μs

In the above code chunk, we are saving the output from running the nichemodel as A_bool and then using the A part of A_bool to construct our initial random network. Within A, 1s indicate an interaction among species and 0s no interaction. In the packages used here, the networks are directed from i (rows) to j (columns), describing the direction of the interaction (i eats j), not of the flow of biomass.

8.1 μs

You can check the connectance of A using:

2.8 μs
co
0.1975
1.2 μs

Here, connectance is calculated as the number of realised links (sum of 1s in A) divided by the number of species in A squared. This end part (species^2) describes the maximum number of possible links in the network A.

7.1 μs

Parameters

Prior to running the BEFW model, you have to create a vector of model parameters using the model_parameters function. Numerous parameter values can be specified within the model_parameters function, however, most of them have default values that are built into the BioEnergeticFoodWebs.jl package. For simplicity, we use the default values here:

7.4 μs
p
Dict
1.0
:e_carnivore
0.85
:Γh
:m_producer
1.0
:c
0.0
:h
1.0
:vertebrates
:tmpA
:rewire_method
:none
:trophic_rank
:w
20×20 Array{Float64,2}:
 0.0   0.0   0.0       0.0       0.0        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   1.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.142857  0.142857  0.142857   …  0.0       0.0       0.0  0.0  0.0
 0.25  0.25  0.25      0.25      0.0           0.0       0.0       0.0  0.0  0.0
 ⋮                                          ⋱  ⋮                             
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.166667  0.166667  0.166667      0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.333333  0.333333  0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0909091     0.0       0.0       0.0  0.0  0.0
:TSR_type
:no_response
:e_herbivore
0.45
:productivity
:species
:efficiency
20×20 Array{Float64,2}:
 0.0   0.0   0.0   0.0   0.0   0.0   0.0   …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.45  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85  …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.45  0.45  0.45  0.85  0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 ⋮                             ⋮           ⋱              ⋮                     
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85  …  0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85     0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.45  0.85  0.85  0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.45  0.0   0.0   0.0  0.0  0.0
:K
1.0
:S
20
:extinctionstime
:is_producer
:dp
#10
:x
:Z
1.0
:extinctions
:dry_mass_293
:np
5
:A
20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0
:ar
:bodymass
:is_herbivore
more
746 μs

For more information and a full list of the parameters and their defaults values type ?model_parameters in the REPL.

If you want to view, check or extract any of the parameter values in p use the p[:name] notation. For example, you can view a vector of each species' trophic rank using:

7.1 μs
439 ns

Simulate

To run the BEFW model, we first assign biomasses at random to each species and then simulate the biomass dynamics forward using the simulate function:

5.0 μs
Dict
:p
Dict
1.0
:e_carnivore
0.85
:Γh
:m_producer
1.0
:c
0.0
:h
1.0
:vertebrates
:tmpA
:rewire_method
:none
:trophic_rank
:w
20×20 Array{Float64,2}:
 0.0   0.0   0.0       0.0       0.0        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   1.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.142857  0.142857  0.142857   …  0.0       0.0       0.0  0.0  0.0
 0.25  0.25  0.25      0.25      0.0           0.0       0.0       0.0  0.0  0.0
 ⋮                                          ⋱  ⋮                             
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.166667  0.166667  0.166667      0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.333333  0.333333  0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0909091     0.0       0.0       0.0  0.0  0.0
:TSR_type
:no_response
:e_herbivore
0.45
:productivity
:species
:efficiency
20×20 Array{Float64,2}:
 0.0   0.0   0.0   0.0   0.0   0.0   0.0   …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.45  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85  …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.45  0.45  0.45  0.85  0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 ⋮                             ⋮           ⋱              ⋮                     
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85  …  0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85     0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.45  0.85  0.85  0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.45  0.0   0.0   0.0  0.0  0.0
:K
1.0
:S
20
:extinctionstime
:is_producer
:dp
#10
:x
:Z
1.0
:extinctions
:dry_mass_293
:np
5
:A
20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0
:ar
:bodymass
:is_herbivore
more
:B
8021×20 Adjoint{Float64,Array{Float64,2}}:
 0.034165   0.175658    0.677504  0.571024   …  0.50232   0.928288   0.0639268
 0.036847   0.0451099   0.434938  0.433516      0.634629  1.12844    0.0804242
 0.0423065  0.0149496   0.318599  0.322718      0.776915  1.32988    0.0978262
 0.0506449  0.00659167  0.26418   0.248923      0.9215    1.4966     0.115019
 0.0617965  0.00361401  0.241252  0.199829      1.05942   1.59266    0.130852
 0.0759696  0.00232627  0.236711  0.165936   …  1.18146   1.60765    0.144321
 0.0935981  0.00168823  0.24454   0.141479      1.27958   1.56103    0.154667
 ⋮                                           ⋱                       
 0.299553   0.200912    0.258188  0.0126103  …  0.0       0.0311834  0.79528
 0.298484   0.200092    0.25727   0.0126249     0.0       0.0311948  0.795905
 0.29731    0.199179    0.256248  0.0126282     0.0       0.0312078  0.796613
 0.296065   0.1982      0.255152  0.0126202     0.0       0.0312223  0.797385
 0.294786   0.197183    0.254016  0.0126011     0.0       0.0312377  0.798198
 0.293516   0.196166    0.252877  0.0125716  …  0.0       0.0312535  0.799025
:t
390 ms

The simulate function requires the model parameters p and the species biomasses bm. In addition, you can specify the timespan of the simulation (using the start and stop arguments), fix a species extinction threshold (using extinction_threshold) and select a solver (using use). For more information type ?simulate in the REPL.

4.9 μs

Output and plot

Once the simulation finishes, the output is stored as a dictionary called out. Within out there are three entries:

  1. out[:p] - lists the parameters

  2. out[:B] - biomass of each species through time

  3. out[:t] - timesteps (these typically increase in 0.25 intervals)

The biomass dynamics of each species can then be plotted. Similar to the DifferentialEquations.jl package, the BioEnergeticFoodWebs.jl package also has it's own built in plotting recipe:

11.3 μs
10.1 ms

You'll notice that the biomass dynamics are noisey during the first few hundred time steps, these are the system's transient dynamics. The dynamics then settle into a steady state where the system can be assumed to be at equilbirum. You'll also notice that some species go extinct and some persist, the initial number of species in the food web (20 in this case) can found using out[:p][:S] and the identity of those that went extinct using out[:p][:extinctions].

The BioEnergeticFoodWebs.jl package also has a range of built in functions that conveniently calculate some of the key metrics of the food web, these include the total biomass, the diversity, the species persistence and the temporal stability:

7.7 μs
biomass
2.398808427844779
62.2 μs
diversity
0.8310935582838079
450 μs
persistence
0.5
249 μs
stability
-0.012570492028439493
166 μs

Each of these functions will output a single value. This value is the average over the last 1000 time steps. For more information, use ? to access the help files on each function in the REPL (e.g., ?species_persistence).

3.5 μs

Variables

Once you've got the BEFW model running, the next step is to vary a variable of interest and rerun. For example, we might be interested in what affect a small change in Z (consumer-resource body mass ratio) has on the estimated food web and its biomass dynamics. The default value for Z is 1.0, but what happens if we increase it to 10.0:

4.5 μs
213 ms

Similarly, what happens if we also increase the carrying capacity (K) of the resource from 1.0 (default) to 5.0:

3.1 μs
555 ms

As you've probably guessed, the main message here is that many variables can be changed in the BEFW model and it's super easy to do so. Some changes will have large effects and some not so much. In the next step, we take this one step further.

6.5 μs

Experiments

The next step is to construct a computional experiment designed to investigate the effect of different variables on population and community dynamics. To do this we construct a gradient of variables as vectors and then simulate the BEFW model multiple times using a loop. To illustrate this, we're going to reproduce example 1 from Delmas et al. 2016. The aim of this example is to investigate the effect of increasing K on food web diversity. In addition, we're also going to allow α (interspecific competition relative to intraspecific competition) to vary and repeat the experiment 5 times with 5 different initial networks.

First, we define the experiment by creating vectors of our variables and fixing the number of repetitions:

5.8 μs
α
645 ns
k
49.6 ms
reps
5
35.0 ns

We then create a dataframe to store the outputs:

2.9 μs
df
αKnetworkdiversitystabilitybiomass
AnyAnyAnyAnyAnyAny
167 ms

and construct a while loop to generate the 5 unique initial networks, each of which contains 20 species with a connectance value of 0.15:

3.2 μs
1.9 ms

We can then run the simulations by looping, using nested for loops, over the unique values of α and K, as well as the 5 unique initial networks. After each simulation we will save each output object to our active project as a JLD2 file and store any output metrics of interest in our dataframe:

3.3 μs
98.6 s

We can then explore the outputs and plot our results. Here, instead of using the built in plotting recipe, we will construct a plot that matches figure 1 in Delmas et al. 2016. Specifically, we will plot food web diversity (y-axis) as a function of K (x-axis) and α (colour):

4.0 μs
variablemeanminmedianmaxnuniquenmissingeltype
SymbolFloat64Float64Float64Float64Int64Int64DataType
1
1.0
0.92
1.0
1.08
3
0
Any
2
:K
2.48181
0.1
1.03291
10.0
10
0
Any
3
:network
3.0
1.0
3.0
5.0
5
0
Any
4
:diversity
0.833252
0.336135
0.843676
1.0
150
0
Any
5
:stability
-0.180435
-1.50381
-0.0165022
-7.14755e-16
150
0
Any
6
:biomass
2.88808
0.301331
2.12444
12.8342
150
0
Any
1.0 s
αKnetworkdiversitystabilitybiomass
AnyAnyAnyAnyAnyAny
1
0.92
0.1
1.0
0.932747
-0.000169469
0.301331
2
0.92
0.16681
1.0
0.938526
-0.000168663
0.426691
3
0.92
0.278256
1.0
0.83412
-0.00014034
0.600914
4
0.92
0.464159
1.0
0.825871
-0.000139997
0.970888
5
0.92
0.774264
1.0
0.82981
-0.27333
1.63807
6
0.92
1.29155
1.0
0.964897
-0.000777308
1.05171
17.0 μs
αKnetworkdiversitystabilitybiomass
AnyAnyAnyAnyAnyAny
1
1.08
0.774264
5.0
0.846921
-0.102008
2.71432
2
1.08
1.29155
5.0
0.880867
-0.358162
3.58376
3
1.08
2.15443
5.0
0.842566
-0.606752
4.05996
4
1.08
3.59381
5.0
0.889399
-0.407231
4.05905
5
1.08
5.99484
5.0
0.649198
-1.16421
5.54389
6
1.08
10.0
5.0
0.999372
-0.00794146
1.8543
14.8 μs
pl
345 ms
shp
900 ns
ls
652 ns
clr
18.2 ms
lbl
867 ns
3.6 s
56.2 ms

Finally, we can save our dataframe as a .csv file:

2.9 μs
"My_data.csv"
682 ms